Semi-Classical Buckling of stiff polymers 
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A quantitative theory of the buckling of a worm like chain based on a semi-classical approximation 
of the partition function is presented. The contribution of thermal fluctuations to the force-extension 
relation that allows to go beyond the classical Euler buckling is derived in the linear and non-linear 
regime as well. It is shown that the thermal fluctuations in the nonlinear buckling regime increase 
the end-to-end distance of the semifiexible rod if it is confined to 2 dimensions as opposed to the 
3 dimensional case. Our approach allows a complete physical understanding of buckling in D — 2 
and in D — 3 below and above the Euler transition. 

PACS numbers: 



I. INTRODUCTION 

During the last few years, the advent of single molecule nanomanipulation [1] has allowed to study the elastic 
properties of DNA and other biopolymers under different physical conditions. In these experiments, the extension of 
single molecule versus an applied stretching force is measured by a variety of techniques including magnetic beads 
[2, 3], optical traps [4, 5], micro-needles [6], hydrodynamic flow [7] and AFM [8]. While the statistical mechanics of 
unconstrained DNA under tension is theoretically well understood in the framework of the Worm Like Chain model 
[9-13] the presence of topological constraints like supercoiling [14-17] or geometrical constraints like protein induced 
kinks and bends [18-22] renders analytical results more difficult. 

Instead of studying the elastic properties of biopolymers under stretching, mechanical properties can also be studied 
by using compression, as long as the chains are smaller than the peersistence length. This has been for example used 
in experiments targeted to measure the force- velocity relation of microtubule growth [23] and in determining the force 
produced by actin filaments [24] . 

With the exception of the work of Odijk [25] who considers a semi-classical evaluation of the partition function in 
the linear regime (Fig. 1 (a)), i.e., below the buckling transition, no calculations have been done that consider the 
non-linear regime of external forces above the critical force (Fig. 1 (b)). Furthermore these calculations are only 
valid well below the transition, although it is the behavior close to the transition on which the force calculations are 
based. In this paper we study thermal fluctuations up to the transition in order to evaluate the scaling of the point 
of buckling with increasing length. 

A computer simulation for 2 and 3 dimensional configurations shows that the thermal fluctuations decrease the 
extension in the buckled state of the polymer in 3 dimensions, but increase it in 2 dimensions (Fig. 4). In this paper 
we show analytically how this happens, by doing a harmonic perturbation calculation around the buckled state. 

As a final note we mention that recently the properties of DNA, like its stiffness and its sequence-specific pairing 
have been exploited to build different kinds of nanostructures [26]. In particular a DNA tetrahedron which has been 
already synthesized could be the building block of extended nanostructures [27]. Our calculations can be used to 
estimate the forces these structures can withstand. 

The paper is organized as follows: we start by describing the geometry and the model used in chapter II, briefly 
treating its classical elastica solutions in chapter III. The main body of the paper consists of a semi-classical calculation 
of the force extension behavior for a WLC with finite length and persistence length below and above the Euler transition 
in chapter IV. We extend these calculations to quartic order below the transition in chapter V, in order to analyze 
the change in buckling transition caused by thermal fluctuations. In chapter VI we compare our calculations with 
simulations. In the concluding chapter VII we discuss our results in the light of several recent experiments with stiff 
biopolymers. 



II. THE PARTITION SUM OF A WORM LIKE CHAIN UNDER COMPRESSION 

We model a stiff polymer as a worm like chain without a twist degree of freedom. In this case, the polymer 
configuration is completely characterized by specifying the unit vector t(s) along the chain, where s is the contour 
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FIG. 1: Force below (a) and above (b) the Euler transition. 



length with < s < L, with L being the chain length. When the chain is submitted to a compressive force T the 
total energy is 



It is custom to write the bending stiffness as A = l P {T)k B T where lp{T) is the orientational persistence length and 
ksT = 1//3 the thermal energy; for example for DNA at room temperature lp w 50nm [28]. All the statistical 
properties of interest can be deduced from the partition function which is a non-trivial quantity to evaluate because 
of the local constraint t 2 (s) = 1 that assures the inextensibility of the chain 



This partition function is nothing but the Euclidean path integral of a quantum particle with mass A, moving on a 
unit sphere under the influence of an external constant force. We are interested in the thermal fluctuations around 
a "classical path". These are easiest to find in polar coordinates. We will fix the force along the a;-axis (see Fig. 
1 (a),(b)) to avoid the chart singularity at the poles. For notational convenience we will choose the polar angle 
t9 € [— 7r/2, 7t/2] such that the uncompressed chain has the coordinates (i?(s), <p(s)) = (0, 0). In these coordinates the 
energy has the form: 




(1) 




(2) 



E[d(s),<p(s)]= f (4 ( cos 2 tf(s)^ 2 (s) +${s)) + T cos tf(s) cos ^(s) 



Jo I 



ds 



(3) 



We can rewrite the energy in dimcnsionlcss variables as: 



<P(t) 




I 



l\ ( cos 2 6(i)c?(i) + e\t) 



+ G 2 



cos#(t) cos cj)(t) > dt 



(4) 



where we have introduced the fluctuation parameter 




and the coupling strength 




The square root in the last expression is in fact the reciprocal of the deflection length [29] of the chain. We are 
interested in the small fluctuation regime and use h as an expansion parameter. The classical path will be the 
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dominant path for h — > 0, i.e., in the rod limit, and thermal fluctuations are taken into account by expanding the 
partition function in fluctuations around this classical path. The partition function will now be a path integral in 
curvilinear coordinates [30]: 



J v 2 [e,<t>]JJW)e- £[e{tum 



(5) 



The determinant of the metric in these coordinates is given by g (9) — cos 2 8. The square root of this determinant, as 
present in the path integral measure, formally takes care of the coordinate independence (chart independence) of the 
measure. It can be understood in a time sliced version, although not without subtleties [30]. This measure term can 
also be formally exponentiated resulting in an extra energy term: 



= J v 2 [0, 



-£[e(t),4>(t)]-e m [6(t),^t)] 



(6) 



with the measure energy term: 



£ m [6(t)} = -6(0) [ dtlogcos 8(t) (7) 
Jo 

The delta function in front of the integral should be understood as being finite using some regularization scheme. The 
classical solutions are obtained through the Euler-Lagrange equations in the next section. We will then proceed by 
incorporating small fluctuations around these classical solutions in section (IV). 

As we will see there are values of the coupling strength where several classical solutions exist with comparable 
Boltzmann weight. These give rise to a bifurcation point in the groundstate. Also since the potential term in (4) is 
positive, there are values of G for which the actual groundstate breaks the rotational symmetry around the direction 
of the applied force. The associated goldstone modes can be excluded by explicitly fixing a direction. 



III. EULER BUCKLING 



In this paper we consider as in Ref. [25] a molecule that has its two ends clamped at fixed orientations 0(0) = 
4>(1) = 8(0) = 8(1) = 0, while the ends can freely move in the plane perpendicular to the force. In the zero fluctuation 
parameter limit the partition function gets only contributions from the classical paths, that minimize the energy. The 
Euler-Lagrange equations are: 

9(t) = -cos8(t)sm8(t)(/) 2 (t) - G 2 sin 8(t) cos 4>(t) (8) 

j (cos 2 8(t)j>(tj) = -G 2 cos8(t)sin(f)(t) (9) 

These equations can be integrated resulting in two classes of solutions: the straight rod solution 

8(t) = <j)(t) = (10) 

and the buckling solutions that read by choosing 8(t) =0 

(t> 2 (t) = 2G 2 (cos 4>(t) - 1 + 2m), to e [0.1) => (11) 
(j)(t) = 2arcsin(Vmsn(£G|m)) => (12) 
cos <f>(t) = 1 - 2m sn 2 (tG\m) (13) 

Here sn() is an elliptic Jacobi function [31]. Solutions with m > 1 are solutions containing loops. They have a higher 
energy in our case. Using the periodicity properties of sn we find for buckling solutions with the boundary condition 
(p(l) = the following relation between m and T 

G = Ly — = 2n K(m) neZ (14) 

Here K(m) is the complete elliptic integral of the first kind. We will label the solutions with n, (f> corresponding to 
the straight rod. Since K(m) is a monotonously increasing function of m we find a smallest force that permits a given 
buckling solution [32]: 

T c = G 2 ^=n 2 ^A = n 2 ^A (15) 
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It is straightforward to calculate the end-to-end distance along the z-axis by integrating the solution along the chain. 
The result is: 

X = L j dt cos (j){t) 

The value of the extension becomes negative under large enough compression. In practice we will be interested in the 
region where the force is small enough that the WLC model is still reasonable. It is easy to see from the buckling 
solution (13) that the compressed chain will be a monotone curve as long asm < 1/2. For higher values of m the 
chain forms an s-shaped curve. 

The energy of the buckling solution is found from (3) to be: 

G 2 ( X 
£ n {m) = — ( 2— + 2m-l 

_ An 2 K 2 (m) / E(m 



.4-^ + 2m -3 (17) 
h \ K(m) / 

with m depending on the force, T, and n through (14). E(m) is the complete elliptic integral of the second kind. 

When comparing the energy of the buckled state with the straight rod configuration we notice that the buckled 
state is always energetically favorable once it is allowed by (15). This transition from straight rod to the buckled 
state is referred to as the Euler transition. When no other constraints are imposed on the solutions the first buckling 
solution, n = 1, will be the favorable solution under compression once the first critical value for the force has been 
reached [32]. When the end of the chain is constrained to be fixed in the origin of the Y"Z-plane, making both ends 
fixed on the z-axis, it is the one loop solution that, when there are no constraints on the rotation of the chain around 
its axes, is the favorable solution. We will for the rest of this article restrict ourself to the unconstrained case. 



IV. SEMICLASSICAL BUCKLING 

For finite values of the fluctuation parameter thermal fluctuations must be taken into account in the evaluation of 
the partition function. We will write the coordinates as: 

6(t) =6 n (t) + 56(t) = S9(t) <j>(t) = 0„(i) + 6<f>(t) (18) 

Here the index n € Z differentiates between the classical solutions (12-14), the straight solution corresponding to 
n = 0. By plugging these relations into the expression for the total energy (including the measure term) we find order 
by order: 

£[e(t),<P(t)]+£ m [6(t)} = 

Jij ^{^" + G 2 COS<? 



— J dt ^j) n 54> — G 2 sin (f> n S(f>} 

~ jf 1 dt |i(<50) 2 - l -G 2 cosUHf + \m 2 \{G 2 cos^ + 0J(50) 2 } 



h . 
+ 1 

+ ■■■ (19) 



The first term is just the energy as given by (17) for the buckle solutions. The second term is zero when we look at 
chains with fixed boundary conditions (Dirichlct boundary conditions). The third term represents the lowest order 
that accounts for thermal fluctuations and is in the focus of our attention. Note that the measure term will only show 
up in the quartic order fluctuations (of order h), since for the Gaussian distribution the fluctuations are of order \fh 
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A. Harmonic fluctuations below the Euler transition 



We first consider the regime below the critical force T c , G < G c = n, where the classical solution is the straight 
rod. The partition function to lowest order around this ground state has the simple form: 

Z = cM-G 2 /h) J V[S6, <50] cxp j-i J 1 dt Q(<56>) 2 - \&{&e?) J cxp j-i £ dt Q(«50) 2 - l -G 2 (<50) 2 ) J (20) 

The resulting path integral is the product of the partition sums, in Euclidean time, of 2 independent harmonic 
oscillators with a frequency squared of — G 2 . When we consider first the azimuth, 0, contribution it is in fact a 
harmonic oscillator on the circle where angles that differ a full period are equivalent. The pathintcgral in that case 
can be expressed as a sum over the harmonic oscillator on the real line by summing over all equivalent end points 
(see e.g. [30] chapter 6): 

^circle(0(O) = 0, 0(1) = 0) = ZlineMO) = 0, 0(1) = 27171) (21) 



n— — oo 



[ G _tf 3 (0,6-2^0 cot G/h) (22) 




V2ttK V sin G 

The elliptic theta function, $3(0,(7) [33] diverges for G = n/2, only half the critical force, which seems to be odd at 
first sight. The reason behind this is that for G — ir/2 all equivalent paths have the same weight, there is no cost in 
increasing the winding number. As a first correction we note that higher order corrections considerably temper the 
potential abyss for larger fluctuations in which case we can neglect the contributions from the winding by taking the 
domain of to be the real line. This results in an improved estimate for the partition sum: 

(23) 



The same kind of reasoning holds for the polar angle. Here we do not have winding, but formally an oscillator 
in a box. Since we again assume the fluctuations to be small it is possible to extend the domain to the real axis. 
Although we have the equivalence (9, 0) ~ (n — 8, + tt) again the results do not hold for larger fluctuations that 
have a weight that does differ substantially from zero. So by taking the polar angle also covering the real line we are 
only overcounting configurations that do not contribute to the path integral. The final result is then: 

* = exp(- G >)^^ (24) 

This partition sum diverges at the caustics, G — tt. Note that this is exactly the critical point for Euler buckling. 
Here it is caused by the harmonic potential being just strong enough to cancel the kinetic term (i.e. the bending 
energy), making large fluctuations favorable and thus invalidating the harmonic approximation. Unlike the $3- 
function divergence here we can not just dismiss these larger fluctuations, since they do not come from a topological 
disconnected region in configuration space and as such are indeed an indication that the groundstate is suffering from 
an instability. 

The force extension behavior is readily obtained, as an approximation, from the partition function: 

= i(l-2^(l-GcotG)) (25) 

This expression diverges again at the Euler transition. Since we have approximated the total extension X = 
L J dt cos 8 cos to quadratic order in the fluctuations around the classical solution, the deviation of above expression 
for the extension from the straight rod actually gives, up to a factor L, the variance of the fluctuations averaged over 
the chain. When this variance is large not only the harmonic approximation to the partition sum breaks but the force 
extension approximation breaks down as well. From these considerations we expect the above force-extension relation 
to hold as long as n — G S> h/2n. From this observation one is tempted to conclude that the rod will start to buckle 
at a force shifted downwards from the Euler transition force following a scaling law for small h of : 



F c ~ Ti 0) (l-Ch) (26) 
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with C a constant of order 1. This is a well known result from Ref. [25]. We will have to adjust this picture when 
taking higher order terms into account, as we will see in section V, because the linear scaling tells us only something 
about the validity of the quadratic approximation. 

For small forces, G -C 1, wc find from (25) for the extension of the chain: 

W=i(l-^(1 + ^)) (27) 
For G = this is the extension of the chain shortened by thermal fluctuations alone. 



B. Harmonic fluctuations above the Euler transition 



The harmonic correction to the classical solution has again the form of an harmonic oscillator, but now with a 
"time" dependent oscillator frequency. The azimuth and polar part of the fluctuation factor again decouple: 

Z = e-^WF+Fo (28) 

The classical solution is given by (17). In principle we should sum over all classical buckling solutions that are allowed 
at a given force. The energy difference is nonetheless big enough that we can neglect the contribution of higher 
buckled configurations. 

The azimuth contribution has the form (after partial integration): 

F = Jv[S4>]exp (-^flj dt5<f>f^ (29) 
with the harmonic fluctuation operator given by 

d 2 

T 4> = --^l -GPcoBfa® (3°) 

where <pi{t) is the classical n = 1 buckling solution. The fluctuation factor can be written using Gaussian integration 
in terms of functional determinants as: 



F rh = 



det(T ) 



d 2 - 1/2 



(31) 



The determinant of the fluctuation operator can be calculated using the Gelfand-Yaglom formula [30] . To do so we 
have to find a solution D${t) of the differential equation 

T^D^(t) = (32) 

with boundary conditions Dj,(Q) = and D<f,(0) = 1. The determinant det(T^) is then given by D^l). Changing 
variables to x = Gt the differential equation has the form of a Lame equation [34] (the Laplacian in ellipsoidal 
coordinates): 



d 2 y(x) 
dx 2 



+ {1 - 2msn 2 (x\m)}y(x) = (33) 



With the given coefficients there exists one double periodic solution (also called Lame polynomial) given by a Jacobi 
elliptic function 

y(t) = cn(Gt) (34) 

This solution has not the right boundary conditions, but using D'Alemberts construction [30], that gives another 
independent solution, we can construct the solution with the right boundary conditions: 

/"* dt' 
D <p (t)=y(t)y(0) 1 



o y 2 (t') 

sn(Gt\m) dn(Gi|m) - E(Gt\m) cn(Gt\m) 
G(l -m) 



+ tcn(Gt\m) (35) 
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Here we adhere to the notation for the Elliptic Integral of the second kind as used in Abramowitz [31], see also 
appendix (A). The function dn is the last Jacobi elliptic function we need. 

With this solution to the Lame equation we find the fluctuation determinant as: 

det (f 



D <t> —rf- = ^(1) 



d 2 , 
dct( -^ 



G(l-m) 

Now we can make use of the relation G = 2K(m) (14) to simplify this result to 



sn(G|m) dn(G|m) - E(G|m ) cn(Gjm) + 



_ E(m)-(l-m)K(m) 
U ^ (l-m)K(m) [6() 



from which we obtain the fluctuation factor: 



(1 - m)K(m) (38) 



27r/i(E(m) - (l-ra)K(ra)) 



Since for small to, E(m) — (1 — to) K(m) ~ to,7t/4, the fluctuation factor diverges at the Eulcr transition. This is not too 
surprising since in the 2 dimensional configuration, there are with forces close to the buckling transition three classical 
solutions with comparable energies with only small barriers in between, allowing larger thermal fluctuations than 
admissible for a harmonic approximation. Would we forbid out-of-plane fluctuations the picture is that fluctuations 
would grow with increasing force just below the Euler transition. Just above the Euler transition the chain will 
fluctuate between the two possible buckled configurations, analogous to quantum tunneling. Finally the buckling will 
stabilize with increasing force to one of the two configurations. 

We now come to the out-of-plane fluctuations. The fluctuation determinants can again be calculated using the 
Gclfand-Yaglom method. We are now looking for a solution of (with x = Gt): 



d 2 y(x) 
dx 2 



+ {l + 4m-6msn 2 (x\m)}y(x) =0 (39) 



This happens to be again a Lame equation with the right coefficients to have a double periodic Lame polynomial as 
solution: 

2/o 0*0 = sn(x) dn(x) (40) 
Since 2/0 (0) = we immediately find for the fluctuation determinant: 



det (f fl ) sn(G) dn(G) 



De:=— V; = ^1^ ^ 

and the partition sum diverges. This is caused by the global rotations around the force direction connecting a 
continuum of groundstates. The buckling solution (13) was chosen to be lying in the rry-plane. Since the energy (as 
well as the pathintegral measure) is invariant under rotations around the x-axis we have a continuum of buckling 
solutions. We can make use of this symmetry by integrating only over paths where the angle 9 averages along the 
chain to zero and then integrating separately over the rotation around the x-axis. This can be done in a consistent 
way using the Faddeev-Popov (FP) method [35] developed to fix internal symmetries in quantum field theory. A 
clockwise rotation of the chain by an angle 7 around the x-axis changes the coordinates on the sphere to: 

cos 9 sin (j) — > cos(0 7 ) sin(</> 7 ) = cos 9 sin 4> cos 7 + sin 9 sin 7 

sin 9 — > sin(0 7 ) = — cos 9 sin <p sin 7 + sin 9 cos 7 (42) 

Now we want to fix the average of the 9 angle, 9 := dt9(t), to zero. We define the FP "determinant" through: 



A FP [6,cf>} / . 
Jo 



dj6(9 1 ) = 1 (43) 
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where the argument of the delta function is the average angle of the by 7 rotated chain. Inserting "1" into the 
partition sum (5) results in: 



/* r 

Z = J d-i Jv 2 [9, <f>] ^jA FP S(e 7 )e- £ ^^ 
= 2ttJv 2 {6,cf>} ^(6)A FP 6(8)e- 



£[e(t),<p(t)] 



(44) 



In the last step we have first performed a trivial change of variable of integration and then made use of the invariance 
under rotation of the energy and of the pathintcgral measure. In fact just the invariance of the combination of the 
measure and the Boltzmann factor would have been enough. The FP determinant can be found from the definition 
(43): 



dt sin 4>{t) 



(45) 



Since we are interested in small thermal fluctuations around the classical solution we can assume the fluctuations to 
be such that J dt sin (f>(t) > for all relevant paths. This apparently does not hold anymore close to the bifurcation 
point. Defining Zo to be the partition sum without the FP term, but including the angle fixing delta function, the 
lowest order contribution of the FP term to the partition sum is: 



Z = [ dt (sin 0(f)) Z 
Jo 

= / disin0i(i)Z o 
Jo 

= / yo(Gt)Z =: Z FP Z 

Jo 



(46) 



The last step follows from the definition of (f>i(t) (13). We now fix the global polar angle in the polar fluctuation 
factor: 



F g := 2nZ FP J V[56]5{56) exp J dt5OT 86^ 



(47) 



To see how this procedure formally gets rid of the divergence we note first that the fluctuation operator, as defined 
on the square integrable functions on [0, 1] that are zero on the boundary, is symmetric and so we can find a real 
orthonormal basis {y n } that diagonalizes the operator. Using this basis we write 60(t) = X)^Lo x nVn (Gt). The 
normalized zero mode eigenfunction is given by yo(Gt) — (J dty 2 (Gt))~ 1 / 2 y (Gt) and the eigenvalues are written as 
A„,e.g. A = 0. We now integrate separately over the zero mode: 



Fg 




Vh 



lim . , 
e^o W 2irhDl 



In the last step we regularized the determinant by adding a small linear term: 

f% :=f 9 + ei 



(48) 



(49) 



in effect shifting all eigenvalues A„ by e to the new values = A„ + e. The resulting determinant is then, in first 
order in e, e times the determinant of the reduced operator defined on the orthogonal complement of the zero mode 
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FIG. 2: Relative extension shift from Eqs. (53) with h = L = 1 



eigenvector, since all other linear terms contain the zero mode eigenvalue. The resulting homogeneous differential 
equation has been solved for a similar case in [36]. The somewhat technical calculation is done in an appendix B. 
The resulting determinant is: 



K(m)3m 

Finally the integral over the zero mode squared is given by: 

1/2 



((1 - m) K(m) - (1 - 2m) E(m)) 



dtyUt) 



/o ' v3 

Combining (38), (48), (50) and (51) we find for the partition sum: 



[(1 - m) K(m) - (1 - 2m) E(m)] 



1/2 



_ £i 2mK(m) 



(1 — m) 



/i 3 [E(m) - (1 - m)K(m) 



(50) 



(51) 



(52) 



It is noteworthy that the partition sum does not diverge at the Euler transition, but goes to zero. By approximating 
the Faddeev-Popov determinant by its classical value we are in fact underestimating the amount of configurations 
the closer we come to the bifurcation point. The force extension corrections to the classical force extension curve 
X (T), equation (16), defined as X = X a + X^ + Xg, with the subscript labeling the fluctuation part that causes the 
extension change, are given by: 



Xa, = 



X e = - 



1 dF^ 

1 dF e 

(3F e dT 



= hL 



m(l — m) 



K(m) 



E(m) — (1 + m) K(m) 



-hL 



16K(m)(E(m) - (1 - m)K(m)) \~E(m) 

(E(m) + 3(1 - m) K(m) 
' 16 K 2 (m)(E(m) - (1 - m) K(m)) 



(1— m)K(m) m(l — m)K(m) 



(53) 



These formula are not too illuminating. Plotting the two corrections (Fig. 2) reveals that the corrections to the 
extension caused by thermal fluctuations have an opposite sign. The out-of-plane fluctuations make the chain slightly 
shorter than the classical solution, as is to be expected. The in-plane fluctuations have the opposite effect. This can 
be understood as the extension change by fluctuations in the straight rod direction to be stronger than fluctuations 
away from the rod solution. 

The total extension again diverges when approaching the bifurcation point, both for the X$ and Xq part seperatcly. 
For the azimuth part the reason behind this is the same as in the straight rod case: near the bifurcation point 
fluctuations increase because the two classical solutions, of positive and negative angle, are close to each other and 
as such a quadratic approximation to the force term is not enough. For the polar angle this is not the case since we 
integrated out the fluctuations to equivalent states, but there the FP term (46) is underestimated: as long as the 
deviation of the expectation value of the end point of the chain (proportional to the FP term) from the straight rod 
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is larger than its fluctuations we can expect that our results hold. Close to the bifurcation point however, we are not 
allowed to drop the absolute value sign by going to Eq. (46) and find a lower bound of the FP term in the order of 
the standard deviation of the end point. 

For small m, approaching the bifurcation point, we find from (53): 

Like below buckling the extension diverges because we make an approximation by taking the extension to be 
— jdplogZ. This is not exact when approximating the potential. For the same reasons as below buckling we can 
expect he results not to hold for large relative extension shifts. 



V. QUARTIC ORDER 

Below buckling it is fairly simple to get a good estimate of the force extension curve up to the Euler transition by 
taking higher order fluctuations into account. Since it is the lowest mode that is responsible for the blowing up of 
the partition sum, approaching the transition, we can significantly improve the calculations by including the quartic 
term for this mode. Quartic terms containing other modes hardly improve upon this. In 2 dimensions the corrected 
partition sum is: 

e-° 2 / h /G(tt 2 -G 2 ) f°° dx ,1,2/2 ^ , aG\ 
V sinG J exp(- — (x 2 (tt 2 - G ) + x 4 —) 



2-irh V sinG J_ oc y/^jrh 2ft 

-&/H /G(tt 2 -G 2 ) 7 ^ / 7 4 \ 

eSTK l O ( 55 ) 



y/2nh V sinG 2tVt7K 4 V 2t V 

with 



G V^—& 

7 = — —Jr — ( 56 ) 



4Vh 2Vh 
From which wc find for the force extension relation: 



3G 2 cotG 1 dj 1 dr _ K 3/4 (^ 7 )+K 5/4 (^) N 7 4 (2 d 1 1 dr 



2G(tt 2 -G 2 ) 2 7 dG rdG y K 1/4 (^r) ' 2t 2 \j dG t dG 



(57) 



These solutions can be continued above the transition, but they start to deviate fast from the exact values. The 
more practical use of these calculations is to make an estimate of the forces a rod can endure, before it collapses. 
Assuming ft <C 1, so that we are close to a buckling type of behavior, we can recognize two separate asymptotic 
regions of behavior, depending on the argument of the modified Bessel functions in (57): 

G 2 < GVh we find as asymptotic behavior: 
2\/2 



x~Hi 



r(l/4): 



i.e. X exhibits a finite negative slope. The decrease of the extension with increasing force is substantial. The 
polymer can be considered to buckle. 

• 7T 2 — G 2 ^> G\fh In this region the decrease of the extension is of order ft, the force extension curve being almost 
flat. There is no buckling yet. 

The crossover region and thus the region where the buckling transition is located, is where this argument is of order 
unity. It is of course not possible to pinpoint a precise transition value, but the scaling of the transition shift follows 
from these observations: the force where the instability appears is shifted by thermal fluctuations according to: 

^ c -^(l-GVft) (59) 
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Classic 
Finite h 
Reduced 
Euler 



FIG. 3: Force extension in 2-d for h = 0.01, 0.05 and 0.1. The corresponding reduced transition forces from Eq. (59), with 
C — 1, are shown by the 3 short dashed vertical lines 



With G of order unity. The results for 3 different values of h are drawn in figure 3 together with the corresponding 
reduced transition forces using C = 1. 

We next consider the 3d case. The contribution from the 9 part alone is the same as for the 4> part, which would 
result in a doubling of the difference from the straight rod. But now we also have a term mixing the two lowest modes. 
The fluctuation part of the partition sum is (apart from a constant): 

G(ir 2 -G 2 ) f°° J [°° J ,1,2, 4 G 2 2x 4 G 2 2 2/ 3G 2 tt 2 
Z fl = -A-^ j ^ dx j ^ dz cxp ( -- ( , 2 A 1 + x 4 — + ^ + Z *- + - y) ) (60) 

As first approximation we can use the expectation value of the square of one of the modes, resulting in a modified 7 
for the other mode given by: 



1 I 1 m / a\ (3G 2 -2tt 2 ) 

' 'tt 2 - G 2 + (z 2 ) (61) 



/-2\ _ / * V 3/4l 2hG 2 ;t^/4\ 2hG 2 ) 1 \ ^ ) a 

\ Z I - \ , ( 7r 2_ G 2 ) 2 1 Q2 TT 2 -G 2 



2Vh 

With: 

K3/4( ^27^)+K 5/4 (^ ^) ^2(. 2 -G 2 ) h 

2K l/4( 2hG 2 > 

The resulting extension is then given by: 

X = X 1 + X^ — L (63) 

This approximation slightly overestimates the contribution of the mixing term close to the transition, where the 
behavior is far from Gaussian. A better result can be obtained by treating the mixing term as a perturbation and 
expanding 60. The resulting series expansion one obtains is: 



G(n 2 -G 2 ) ^ 1 f 3G 2 -2tt 2 Y t T(^) <E(^, 1,^) - 2 7 ^ T(^p) j>(^p, f , ^) 
Zfl= sinG 8^J 2r( 2 «+3)/2 ( 64 ) 



with ^(x, y, z) Kummer's function (confluent hypergeometric function). This series converges relatively fast just below 
the Euler transition and one can get a good approximation of the extension below the transition force. For practical 
purposes the first approximation is good enough to characterize the transition shift. It scales with increasing length 
in the same way as the 2d case. We will next compare the predictions of the force-extension realtions, Eqs. (25), 
(53), (58), (63) and (64) with simulations. 
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G 2 

FIG. 4: Comparison of the analytical force-extension with simulations for L — 49 and h = 0.8. The unit of length is the bond 
length in the simulation. 

VI. COMPARISON WITH THE SIMULATION 

The simulations were done with bond and bead models, consisting of 50 beads joined by either a bond consisting 
of a strongly repulsive Lennard- Jones potential and an attractive FENE potential [37, 38], or with a stiff harmonic 
bond. The bending potential was implemented by a cosine angular energy term with a magnitude chosen so that the 
persistence length was comparable to the chain length (L = 49 in terms of the bond length) . The backbone stretching 
parameters were chosen in such a way that possible fluctuations of the bond length can be neglected compared to the 
bending fluctuations. The incxtcndiblc worm-like chain can then be expected to be a reasonable approximation to 
the simulated chain. The general features of the simulation did not depend on the type of bond chosen. 

The simulations and theoretical calculations are plotted in Fig. 4. The value of h was taken rather high in order 
to have a more pronounced fluctuation contribution. The length scale is chosen such that the bond length in the 
simulation model is 1. The 3-d quartic curve was calculated using the modified quartic term. 

The semi-classical results are in good agreement with the simulation data in the region where a semi-classical 
approximation is expected to be valid. It is noteworthy that the increase in extension as predicted by the calculations 
is indeed the same as observed in the simulation. In 2d the quartic corrections below buckling show even for relatively 
large values of h, good agreement with simulations. In 3d, using the simplified approach of modifying the quartic 
interaction to account for the mode mixing (63) the reliability of the calculations close to the Euler transition decreases, 
although the qualitative behavior seems to be good enough for practical purposes. Better results one gets using a 
perturbation expansion (64). The 3-d quartic series curve was calculated using this expansion with the first 20 terms. 
Note though that this last calculation was stopped slightly below the transition force, since it does not converge at 
the transition. 

The effect of the bond length not being fixed is indeed small enough compared to the thermal fluctuations. The 
errorbars are caused by the finite number of simulation rounds. 
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lp 


L 


ft 


A 0) 


T c 


DNA 


50nm 


lOnm 


0.2 


21pN 


llpN 


actin 


9/im 


1.2/im 


0.13 


0.26pN 


0.16pN 


actin with Phalloidin. 


18pim 


0.75/ttm 


0.04 


1.3pN 


l.OpN 


Microtubule 


3.3mm 


9.4^m 


0.028 


1.5pN 


1.4 pN 



TABLE I: reduction of the force needed to buckle for some biopolymers with finite length. The reduction is calculated from 
Eq. (59) with C = 1 

VII. DISCUSSION 

The parameter that determines whether a buckling transition is present is the ratio h of length and persistence 
length of the wormlikc chain. One can roughly say that a buckling transition appears for ratios clearly smaller than 
1. But it is crucial that one takes into account the shift of the apparent transition when a force is extracted from 
the onset of buckling. To illustrate the importance of thermal fluctuations we will discuss the influence they have in 
interpreting data from recent experiments with important biopolymers. Table I shows the persistence length of the 3 
polymers, ds-DNA, actin and Microtubule together with some of the typical lengths and associated transition forces. 
The shifted transition force is calculated from Eq. 59 with (7=1. 

The DNA tetraheda synthesized by Goodman et.al. [27] have sides made of double stranded DNA of a length below 
lOnm. As can be read of from the table, for a lengths of lOnm the force the structure can endure is strongly reduced 
by thermal fluctuations. This has to be taken into account when designing nanostructures based on DNA. 

F-actin is one of the main building blocks of the cytoskeleton. It has a persistence length in the order of 9 — 18/im 
[39] (the higher value is in presence of the toxin Phalloidin). actin can produce forces through polymerization. The 
maximum force it can produce, the stall force, was determined, by Kovar et.al. [24], by measuring the shortest length 
of actin that showed buckling, when growing in between 2 fixed points. The lengths where this was observed are given 
in row 2 and 3 of the table. The force calculation based on classical buckling considerably overestimates the force 
needed to buckle for the measured length since it does not take the thermal fluctuations into account. 

The other important structures in the cytoskeleton are microtubules, hollow highly regular assemblies of filaments, 
having persistence lengths in the order of several mm's [39], the precise value depending on several factors, like the 
growth speed [23] and perhaps the contour length [40]. In buckling experiments by Janson et.al. [23], where the 
growth rate dependence on the applied force was studied, the lengths were such that in this case the shift by thermal 
fluctuations is neglegiblc. 

Nevertheless, the increase of thermal fluctuations when approaching buckling can also be observed in this case. These 
thermal fluctuations increase sharply just before buckling, followed by a strong damping of these fluctuations with 
increasing length (and thus increase of buckling) of the microtubule. Both these effects follow from our calculations. 

The damping of the fluctuations after the onset of buckling can be inferred from the approach of the semiclassical 
solution towards the "zero temperature" classical solution. Below buckling the end point fluctuations increase from 
(y 2 ) = L 2 h/3, the classic result which follows from (27), for a chain with one free end to (y 2 ) w 0.2L 2 \/h for an applied 
force corresponding to the Euler transition, as follows from (58). It should be noted that the geometry of the setup 
in those experiments is not immediately comparable to our calculations since the microtubule in those experiments 
have one end of the chain more or less hinged in a fixed position, the resulting buckling force can be up to a factor 4 
larger than in our case. Qualitatively though the results are comparable and for typical values of a persistence length 
of 3.3mm and a chain length of 20^m we expect the mean fluctuation of the end point to be amplified by a factor 
w 7. This indeed seems to be approximately the case, although a precise analysis of their measurements is outside of 
the scope of this paper. 

Finally, a remarkable result of our calculations is the increase of end-to-end distance by thermal fluctuations of the 
buckled polymers, especially in 2 dimensions. In dense networks of actin filaments confined to the cell cortex, the 
buckling is approximately 2-dimensional. The lengthening of the buckled polymer causes then an apparent stiffening 
of the compressed network by the fluctuations. 

APPENDIX A: ELLIPTIC FUNCTIONS 

The elliptic integrals and the Jacobi elliptic functions are functions of two variables and in the case of elliptic 
integrals of the third kind three. There are different equivalent choices of pairs and the choice generally depends on 
the situation at hand. See also [31, 33]. 



14 



Throughout this paper we use the Jacobi form (with one exception). In that form the variables are called the 
argument, x, and the parameter, m € [0,1]. In the literature the latter is sometimes replaced by the modulus, 
k = \fm. The two variables are separated by a vertical line like in E(x|m). An alternative form is the trigonometric 
form where the variables are the Jacobi amplitude, <f> = am(x|m) and the modulus a defined through sin 2 (a) := m. 
In that case the variables are separated by a backslash. So in the notation that we use we have the elliptic integrals 
of the first, second and third form written as: 

F((f)\m) = F(<t>\a) E(x\m) = E(</>\a) U(n; x\m) = U(n; cj)\a) (Al) 

The integral of the first kind is an exception since it is in fact the inverse of the amplitude function and so F(x\m) 
is identical to x. The complete integral of the first kind is defined as the value of F evaluated at an amplitude of 
7r/2: K(m) := F(7r/2|m). The same holds for the other complete integrals, but now we can make use of the fact that 
am _1 (7r/2|m) = K(m) and so: 

E(m) := E(K(m)|m) U(n\m) := U(n;K(m)\m) (A2) 

The double periodic Jacobian elliptic functions are defined as: 

sn(x\m) := sin(am(x|m)) cn(x|m) := cos(am(x|m)) dn(x\m) := -^-am(x|m) (A3) 



APPENDIX B: GENERALIZED LAME EQUATION[36] 

We are looking for a solution of the generalized Lame equation: 

y+p(x)y = (Bl) 

where p(x) = 1 + Am — e — 6msn 2 (x) and we are especially interested in the small e limit. The product M(x) = 
yi(x)y2(x) of two solutions satisfies the third order differential equation: 

M + ApM + 2pM = (B2) 

We will now construct a solution of this last equation as a series in sn(x). Write M = J2 ns o a ™ snTl ( x )- Substitution 
leads to the following relation between the coefficients: 

a n m(n 3 + 3n 2 - 22n - 24) + o n+2 (4(l + Am - e)(n + 2) - (n + 2) 3 (1 + m)) + a n+4 (n + 4)(n + 3)(n + 2) = (B3) 

To get a finite number of terms, the highest power has to be 4 and we find as solution: 

M(x) = 9m 2 sn 4 (x) - 3m(3 + e) sn 2 (x) + 3e(l - m) + e 2 (B4) 

Suppose yi(x) is one of the 2 solutions of (Bl) that make up M. We can use the D'Alembert construction to get 
another independent solution so that y 2 can be written as (the Wronskian is constant): 

y 2 (x) = B Vl (x) + C yi (x) f dx'^— (B5) 

Jo Vi\ x ) 

Using the definition of M{x), and assuming M(x) to be positive, we can express y(x) in terms of M(x) as: 

Inserting this function into the Lame equation results in: 

2M(x)M{x) - M 2 {x) + C 2 + 4(1 + Am - e - 6m sn 2 (x\m)) M 2 (x) = (B7) 
and the C(ra) is found by inserting the solution for M(x) (B4): 

C(m) = ±2^e(3m(3 + e) - e(l + Am - e)(3(l - m) + e))(3(l - m) + e) = ±6^\/3m(l - m) (B8) 
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The integrand in the exponential of (B6) has poles at 

_ 3 + e ± a/9 - 6e(l - 2m) - 3e 2 ^ 3 + e ± (3 - e(l - 2m)) 
± 6m 6m 

Since sn 2 (x) e [0, 1], by choosing e < we can force the integrand to be regular for the parameter m e [0, 1]. C(m) 
on the other hand is now imaginary , so that we have to look at linear combinations of the two solutions for a real 
valued solution of the homogeneous equation. Noting that M(x) is now strictly negative we find as the solution with 
the proper boundary conditions: 

2J-M(fi) , — . . f-\C{m)\ f x dx' \ 
We will now evaluate the integral appearing in the sinus: 



/' 

Jo 



dx' -l r dx' i r dx' 

+ 



o M(x') 9m 2 (p + — p~)p+ J 1 — l/p+ sn 2 (x'|m) 9m 2 (p + — p~)p~ J 1 — l/p~ sn 2 (x'|m) 

= Qm 2 ^-^)^ ^ 1 ^^'-^ ^^-^)^ ^ 1 ^^'^ ^ 

The elliptic integral of the third kind, II(n; x\m), has a behavior that depends on the value of the characteristic n [31]. 
In our case we have characteristics l/p+ = m(l — em/3) € (m, 1) and l/p~ = 3m/(e(l — m)) < both corresponding 
to so called circular cases. 

The functional determinant is then given by: 



^-n(l/p+|m) - ^n(l/p_|m) 



(Bll) 



K(m)\C(m)\ \9m 2 ( P+ -p_) 
where we introduced the complete elliptic integral of the third kind, defined as usual from the incomplete one: 

Yl(n\m) := IT(n; K(m)|ra) = hl(n; 2K(m)|m) (B12) 

The integral with the l/p+ characteristic can be easily evaluated by expanding around characteristic m, leading to: 

/■ K(m) dx 1 

U(l/ P+ \m) = / + O(c) - E(m) (B13) 

Jo dn (x\m) 1-m 

Note that we only need up to first order in e and the factor multiplying the sinus in Eq. (Bll) goes as M(0) ~ e. 

The integral with the negative characteristic can be written as an elliptic integral with positive characteristic using 
the following definition [31]: 

1_^ (l-m) 2 

\-p- 3m V ' 

resulting in: 

n(l/p_|m) - -p-(l-m) n(jv_ Irre) - P ' TO K(m) 
v /F 1 y (l-p_)(l-mp_) v 1 y l-p_m v ; 

= -p-(l - m)II(A^_|m) -p_mK(m) (B15) 

The characteristic is less trivial, since we can't expand around characteristic 1, where the elliptic integral 
diverges. We can express the integral in terms of yet another elliptic function, Hcuman's Lambda function A (z|m), 
as: [31] 

U(N- \m) = K(m) + ^n5[l - A (z\m)} (B16) 

where 



III 3m 



(l-N)(N-m) V^VVl 1 -™) 3 



The Lambda function we can express, again following Ref. [31] in other elliptic functions as: 



U(N.\m) = KM - _ E(m) + ^ _ — + 0(V=i) (B19) 



And finally we find for the functional determinant: 



K(m)3m 
which is Eq. (50) 
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A (z\m) = - (K(m) E(z|l - m) - (K(m) - E(m)) F(am(z)|l - m)) (B18) 

7T 

Using the small argument expansion of the elliptic integrals: E(z|l — m) = z wc find : 



D « " -^Th^rV sin f a 'n^ 1 \ [K(m)(l - m) - (1 - 2m) E(m)] + tt 
K(to)|6(to)| \9m(l — m) 

[K(m)(l-m)-(l-2m)E(m)] (B20) 
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